An Approach to Modeling Coupled Thermal-hydraulic-chemical Processes in Geothermal Systems

Interactions between hydrothermal fluids and rock alter mineralogy, leading to the formation of secondary minerals and potentially significant physical and chemical property changes. Reactive transport simulations are essential for evaluating the coupled processes controlling the geochemical, thermal and hydrological evolution of geothermal systems. The objective of this preliminary investigation is to successfully replicate observations from a series of hydrothermal laboratory experiments [Morrow et al., 2001] using the code TOUGHREACT. The laboratory experiments carried out by Morrow et al. [2001] measure permeability reduction in fractured and intact Westerly granite due to high-temperature fluid flow through core samples. Initial permeability and temperature values used in our simulations reflect these experimental conditions and range from 6.13 × 10 to 1.5 × 10 m and 150 to 300 °C, respectively. The primary mineralogy of the model rock is plagioclase (40 vol.%), K-feldspar (20 vol.%), quartz (30 vol.%), and biotite (10 vol.%). The simulations are constrained by the requirement that permeability, relative mineral abundances, and fluid chemistry agree with experimental observations. In the models, the granite core samples are represented as one-dimensional reaction domains. We find that the mineral abundances, solute concentrations, and permeability evolutions predicted by the models are consistent with those observed in the experiments carried out by Morrow et al. [2001] only if the mineral reactive surface areas decrease with increasing clay mineral abundance. This modeling approach suggests the importance of explicitly incorporating changing mineral surface areas into reactive transport models. INTRODUCTION The search for sources of renewable energy has promoted interest in geothermal energy exploration and development. Interactions between hydrothermal fluids and rock alter the primary mineralogy, leading to formation of secondary minerals and potentially significant physical and chemical property changes [Xu and Pruess, 2001; Xu et al., 2003]. Reactive transport simulations provide an objective and systematic means of investigating and evaluating the competitive interactions between fluids and rockforming minerals over space and time for a range of conditions [Steefel et al., 2005]. However, endeavors to model coupled thermal-hydrologic-chemical (THC) processes are complicated by the number of poorly constrained processes that occur under hydrothermal conditions. In order to maximize their utility, reactive transport models must adequately integrate the various fundamental processes that control system properties and evolution. A major challenge in the field of reactive transport modeling has been the development of general models that will transfer to different geothermal sites and conditions. The many assumptions, uncertainties, and approximations intrinsic to reactive transport models restrict our ability to draw general conclusions regarding the behavior of natural geothermal systems. Although hydrothermal fluids and rocks interact under a diverse range of conditions, most studies necessarily present simulation results that are specific to a particular set of conditions and parameters [e.g., Xu and Pruess, 2001; Xu et al., 2001; Maher et al., 2006]. If reactive transport models are to be more broadly applied to problems such as the analysis and exploration of geothermal systems, then these models must be transferable from one system to another (i.e., not site specific). This requires modeling techniques that are sufficient enough to apply over a wide-range of conditions. Reactive transport models rely on mathematical formulations to describe the coupled THC processes occurring in hydrothermal environments [Johnson et al., 1998]. The general mathematical expression typically used to predict mineral dissolution or precipitation rates is based on transition state theory and has the form [Aagaard and Helgeson, 1982; Helgeson et al., 1984; Lasaga, 1984]: rate = A⋅ k⋅ f ∆G ( ) (1) where A is the specific reactive surface area, k the rate constant, and f (∆G) a function which gives the dependence of the rate on the Gibbs free energy (∆G). Equation (1) illustrates the overall rate dependence on: (i) the proximity to equilibrium (i.e., the degree of saturation, which is defined in terms of the Gibbs free energy), (ii) rate constants, and (iii) mineral surface areas. Changes to any of these terms can contribute to significant differences in the predicted mineral dissolution/precipitation behavior and, by extension, the evolution and behavior of the system. Comparisons of various studies [e.g., Martin and Lowell, 1997; Johnson et al., 1998; White and Brantley, 2003; André et al., 2006; Maher et al., 2006; Yasuhara and Elsworth, 2006] indeed indicate that differences in kinetics, surface area, solution chemistry, and mineral composition and solubility can cause significant differences in predicted mineral precipitation and permeability evolution. Therefore, in addition to appropriate mathematical formulations of the coupled processes, reactive transport models require accurate compositional, thermodynamic, and kinetic data. Unfortunately, values for parameters such as mineral compositions and solubilities, fluid chemistry, rate constants, and reactive surface areas are often difficult to constrain or unavailable [White and Peterson, 1990; Alekseyev et al., 1997; Cama et al., 2000; Brantley, 2003; Maher et al., 2006]. The multitude of unconstrained variables presents a significant challenge. For the purposes of this study, we have reduced the problem to a more tractable size by focusing on the numerical treatment of reactive surface areas. Although all the parameters mentioned above introduce uncertainty into reactive transport models, mineral surface areas remain one of the most uncertain and poorly quantified parameters [White and Brantley, 1995; Bethke, 1996; Alekseyev et al., 1997; Cama et al., 2000; Steefel, 2001; Lasaga and Lü ttge, 2003; Xu et al., 2004a; Maher et al., 2006; Brantley, 2003, 2008]. Relatively few studies have attempted to simulate physical changes in mineral surface areas or to evaluate their effects on dissolution and precipitation rates [Sonnenthal and Ortoleva, 1994; Brantley and Conrad, 2008]. Surface areas are notoriously difficult to ascertain [Nagy et al., 1991; Nagy and Lasaga, 1992; Stillings and Brantley, 1995; Amrhein and Suarez, 1992; Ganor et al., 1995, 1999; Cama et al., 2000] and are known to change with time due to poorly quantifiable phenomena. Processes affecting reactive surface area include coating of reactive mineral phases by other, less reactive minerals (armoring), dissolution pitting and etching, and creation of isolated porosity through crack healing and sealing [Brantley et al., 1990; Chester et al., 1993, Brantley, 2003; Brantley and Conrad, 2008]. The findings of White and Brantley [2003] and Maher et al. [2006] indicate that these changes in mineral surfaces may significantly contribute to the observed changes in rates. Using the popular fluid flow and geochemical transport code TOUGHREACT [Xu and Pruess, 1998; Xu et al., 2004a], we present a relatively straightforward method for representing the progressive loss of reactive surface area due to occlusion by secondary precipitates. Since a critical benchmark for any numerical model is accurate simulation of well-constrained experiments, model success is predicated on agreement with a series of flow-through laboratory experiments conducted on cylindrical samples of Westerly granite at temperatures from 150 to 300 °C [Morrow et al., 2001]. Both intact and fractured granitic samples were used in the experiments, resulting in a range of initial permeability values. The objective of this preliminary investigation is to evaluate the use of a one-dimensional model that explicitly considers the temporal evolution of reactive surface areas for simulating the physical and chemical evolution of fractured granite under various hydrothermal conditions. The results presented in this paper can be compared to results from more complex multidimensional models. The intent is to better quantify and characterize the processes dominating the evolution of natural geothermal systems. Because rates are controlled by a combination of factors, an improved mechanistic or process-based understanding should facilitate the future development and utility of these models for field applications. SURFACE AREA EFFECTS Modeling reaction rates and their subsequent effects on system evolution requires knowledge of the total mineral surface areas in contact with the aqueous phase. The surface area, A, used in Eq. (1) typically is based on either gas adsorption measurements (BET) or geometric estimates [Maher et al., 2006; Brantley and Conrad, 2008]. However, BET surface areas are consistently higher than those based on geometric estimates [Dorn, 1995; Brantley et al., 1999; White and Brantley, 2003]. Additionally, there are substantial uncertainties associated with both Table 1: Summary of experiments from Morrow et